#include "stdafx.h"
#include "matrix.h"
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <iterator>
#include <algorithm>
#include <cmath>

template <typename T> 
Matrix<T>::Matrix(int N)
{
	itemsCount = N*N;
    items = items = std::shared_ptr<T>(new T[itemsCount], []( T *p ) { delete[] p; });
	rows = columns = N;
    rowOffset = 0;
    columnOffset = 0;
    parent = NULL;
};

template <typename T> 
Matrix<T>::Matrix(int N, int M)
{
	itemsCount = N*M;
    items = std::shared_ptr<T>(new T[N*M], []( T *p ) { delete[] p; });
	rows = N;
	columns = M;
    rowOffset = 0;
    columnOffset = 0;
    parent = NULL;
};

template <typename T> 
Matrix<T>::Matrix(Matrix<T> &parent, int N, int M)
{
    itemsCount = N*M;
    items = parent.items;
	rows = N;
	columns = M;
    this->rowOffset = 0;
    this->columnOffset = 0;
    this->parent = &parent;
}

template <typename T> 
Matrix<T>::Matrix(Matrix<T> &parent, int rowOffset, int columnOffset, int N, int M)
{
    itemsCount = N*M;
    items = parent.items;
	rows = N;
	columns = M;
    this->rowOffset = rowOffset;
    this->columnOffset = columnOffset;
    this->parent = &parent;
}

template <typename T> 
Matrix<T>::~Matrix(void)
{
};



template <typename T> 
void Matrix<T>::Initialize(T number)
{
	for(Matrix<T>::iterator i = this->begin(); i != this->end(); i++)
		{
			(*i) = number;
		}
};

template <typename T> 
void Matrix<T>::Initialize(int min, int max)
{
	for(Matrix<T>::iterator i = this->begin(); i != this->end(); i++)
	{
		(*i) = rand() % max + min;
	}
};

template <typename T> 
void Matrix<T>::Initialize(T* collection)
{
	items = std::shared_ptr<T>(collection, []( T *p ) { delete[] p; });
};

template <typename T> 
void Matrix<T>::Zeros()
{
	Initialize(0);
};

template <typename T> 
void Matrix<T>::Ones()
{
	Initialize(1);
};

template <typename T> 
void Matrix<T>::Eyes()
{
	if(columns != rows)
		throw INVALIDSIZE;

	Initialize(0);
	int a = 0;
	for(Matrix<T>::iterator i = this->begin(); i != this->end(); i++)
	{
		if(a++%(columns+1) == 0)
		{
			(*i) = 1;
		}
	}
};

template <typename T>
template <typename T2>
Matrix<T2> Matrix<T>::Convert() const
{
    Matrix<T2> copy(rows, columns);

    int index;
    int copyIndex;
	for(int n=0;n<rows;n++)
	{
        for(int m=0;m<columns;m++)
        {
            index = CalculateIndex(n, m);
            copyIndex = copy.CalculateIndex(n, m);
		    copy.items.get()[copyIndex] = (T2)items.get()[index];
        }
    }

    return copy;
}

template <typename T> 
Matrix<T> Matrix<T>::Clone() const
{
    Matrix<T> &copy = this->Convert<T>();
    return copy;
};

template <typename T> 
int Matrix<T>::CalculateIndex(int N, int M) const
{
    if (parent == NULL)
        return N * columns + M;

    return (N + rowOffset) * parent->columns + M + columnOffset;
};

template <typename T> 
T Matrix<T>::GetMinorForElement(int N, int M) const
{
	int i=0;
	int row=0;
	int column=0;
	T value;
	T* values = new T[(rows-1)*(columns-1)];
	
	for(row=0;row<rows;row++)
	{
		if(row == N)
			continue;

		for(column=0;column<columns;column++)
		{
			if(column == M)
				continue;

			value = ElementAt(row, column);
			values[i] = value;
			i++;
		}
	}

	Matrix<T> matrix(columns-1, rows-1);
	matrix.Initialize(values);

	T det = matrix.CalculateDeterminant();

	return det;
}

template <typename T> 
Matrix<T> Matrix<T>::GetComplementsMatrix() const
{
	Matrix<T> complements(columns, rows);

	int converter = 0;
	T minor = 0;
	
	for(int row=0;row<rows;row++)
	{
		for(int column=0;column<columns;column++)
		{
			converter = pow((double)-1, (double)(column+row));
			minor = GetMinorForElement(row, column);
			complements.SetElementAt(row, column, converter * minor);
		}
	}

    return complements;
}

template <typename T>
void Matrix<T>::Show(void)
{
	cout<<(*this);
};

template <typename T>
int Matrix<T>::GetRowsCount() const
{
    return rows;
}

template <typename T>
int Matrix<T>::GetColumnsCount() const
{
    return columns;
}

template <typename T> 
T Matrix<T>::ElementAt(int N, int M) const
{
	int index = CalculateIndex(N, M);
	return items.get()[index];
};

template <typename T> 
void Matrix<T>::SetElementAt(int N, int M, T value)
{
	int index = CalculateIndex(N, M);
	items.get()[index] = value;
};

template <typename T>
Matrix<T> Matrix<T>::Transpose() const
{
    Matrix<T> transposed(columns, rows);

	for(int n=0;n<rows;n++)
	{
        for(int m=0;m<columns;m++)
        {
			transposed.SetElementAt(m, n, ElementAt(n, m));
        }
    }

    return transposed;
}

template <typename T>
Matrix<double> Matrix<T>::Inverse() const
{
	T determinant = CalculateDeterminant();

	if(determinant == 0)
		throw INVALIDOPERATION;

	Matrix<T> complements = GetComplementsMatrix();
	Matrix<T> transposed = complements.Transpose();
	Matrix<double> transposedDouble = transposed.Convert<double>();
	
    return transposedDouble * (1.0/(determinant*1.0));
}

template <typename T> 
T Matrix<T>::CalculateDeterminant() const
{
	if(columns != rows)
		throw INVALIDSIZE;

	if(columns > 3)
		throw UNSUPPORTEDSIZE;
	else if(columns == 1)
		return ElementAt(0,0);
	else if(columns == 2)
		return ElementAt(0,0) * ElementAt(1,1) - ElementAt(0,1) * ElementAt(1,0);
	else 
	{
		T part1 = ElementAt(0,0) * ElementAt(1,1) * ElementAt(2,2),
			part2 = ElementAt(0,1) * ElementAt(1,2) * ElementAt(2,0),
			part3 = ElementAt(0,2) * ElementAt(1,0) * ElementAt(2,1),
			part4 = ElementAt(0,2) * ElementAt(1,1) * ElementAt(2,0),
			part5 = ElementAt(0,1) * ElementAt(1,0) * ElementAt(2,2),
			part6 = ElementAt(0,0) * ElementAt(1,2) * ElementAt(2,1);

		T det = (part1 + part2 + part3) - (part4 + part5 + part6);
		return det;
	}

	return 0;
}

template <typename T>
T Matrix<T>::Min()
{
	Matrix<T>::iterator last = this->end();
	last -=1;
	Matrix<T>::iterator it = std::min_element(this->begin(), last);
	return (*it);
}

template <typename T>
T Matrix<T>::Max()
{
	Matrix<T>::iterator last = this->end();
	last -=1;
	Matrix<T>::iterator it = std::max_element(this->begin(),last);
	return (*it);
}


template <typename T>
Matrix<T>& Matrix<T>::operator=(Matrix<T>& otherMatrix)
{
	if(columns != otherMatrix.columns || rows != otherMatrix.rows)
		throw INVALIDSIZE;

	Matrix<T>::iterator j = otherMatrix.begin();

	for(Matrix<T>::iterator i = this->begin(); i != this->end(); i++)
	{
		(*i) = (*j);
		j++;
	}
	
	return *this;
};

template <typename T>
bool Matrix<T>::operator==(Matrix<T>& otherMatrix)
{
	if(columns != otherMatrix.columns || rows != otherMatrix.rows)
		return false;

	Matrix<T>::iterator j = otherMatrix.begin();

	for(Matrix<T>::iterator i = this->begin(); i != this->end(); i++)
	{
		if((*i) != (*j))
			return false;

		j++;
	}
	
	return true;
};



template <typename T>
Matrix<T> Matrix<T>::operator+(Matrix<T>& otherMatrix)
{
	if(columns != otherMatrix.columns || rows != otherMatrix.rows)
		throw INVALIDSIZE;

	Matrix<T> result(rows, columns);

	Matrix<T>::iterator i =  this->begin();
	Matrix<T>::iterator j = otherMatrix.begin();

	for(Matrix<T>::iterator k = result.begin(); k != result.end(); k++)
	{
		(*k) = (*i) + (*j);
		i++;
		j++;
	}
	
	return result;
};

template <typename T>
Matrix<T> Matrix<T>::operator-(Matrix<T>& otherMatrix)
{
	if(columns != otherMatrix.columns || rows != otherMatrix.rows)
		throw INVALIDSIZE;

	Matrix<T> result(rows, columns);

	Matrix<T>::iterator i =  this->begin();
	Matrix<T>::iterator j = otherMatrix.begin();

	for(Matrix<T>::iterator k = result.begin(); k != result.end(); k++)
	{
		(*k) = (*i) - (*j);
		i++;
		j++;
	}
	
	return result;
};

template <typename T>
Matrix<T> Matrix<T>::operator*(Matrix<T>& otherMatrix)
{
	return MatrixMultiplication(*this, otherMatrix);
};

template <typename T>
Matrix<T> Matrix<T>::operator/(Matrix<T>& otherMatrix)
{
	if(columns != otherMatrix.columns || rows != otherMatrix.rows)
		throw INVALIDSIZE;

	throw NOTIMPLEMENTED;

	Matrix<T> result(rows, columns);

	//todo:
	
	return result;
};



template <typename T>
Matrix<T>& Matrix<T>::operator+=(Matrix<T>& otherMatrix)
{
	if(columns != otherMatrix.columns || rows != otherMatrix.rows)
		throw INVALIDSIZE;

	Matrix<T>::iterator j = otherMatrix.begin();

	for(Matrix<T>::iterator i = this->begin(); i != this->end(); i++)
	{
		(*i) += (*j);
		j++;
	}
	
	return *this;
};

template <typename T>
Matrix<T>& Matrix<T>::operator-=(Matrix<T>& otherMatrix)
{
	if(columns != otherMatrix.columns || rows != otherMatrix.rows)
		throw INVALIDSIZE;

	Matrix<T>::iterator j = otherMatrix.begin();

	for(Matrix<T>::iterator i = this->begin(); i != this->end(); i++)
	{
		(*i) -= (*j);
		j++;
	}
	
	return *this;
};

template <typename T>
Matrix<T>& Matrix<T>::operator*=(Matrix<T>& otherMatrix)
{
	if(this->parent != NULL)
		throw CHILDMATRIX;
	
	*this = MatrixMultiplication(*this, otherMatrix);

	return *this;
};

template <typename T>
Matrix<T>& Matrix<T>::operator/=(Matrix<T>& otherMatrix)
{	
	throw NOTIMPLEMENTED;

	//todo:

	return *this;
};



template <typename T>
Matrix<T>& operator+=(Matrix<T>& matrix, T value)
{
	for(Matrix<T>::iterator i = matrix.begin(); i != matrix.end(); i++)
	{
		(*i) += value;
	}
	
    return matrix;
};

template <typename T>
Matrix<T>& operator-=(Matrix<T>& matrix, T value)
{
	for(Matrix<T>::iterator i = matrix.begin(); i != matrix.end(); i++)
	{
		(*i) -= value;
	}
	
    return matrix;
};

template <typename T>
Matrix<T>& operator*=(Matrix<T>& matrix, T value)
{
	for(Matrix<T>::iterator i = matrix.begin(); i != matrix.end(); i++)
	{
		(*i) *= value;
	}
	
    return matrix;
};

template <typename T>
Matrix<T>& operator/=(Matrix<T>& matrix, T value)
{
	for(Matrix<T>::iterator i = matrix.begin(); i != matrix.end(); i++)
	{
		(*i) /= value;
	}
	
    return matrix;
};



template <typename T>
Matrix<T> operator+(Matrix<T>& matrix, T value)
{
    Matrix<T>& copy = matrix.Clone();
    copy += value;

    return copy;
}

template <typename T> 
Matrix<T> operator-(Matrix<T>& matrix, T value)
{
    Matrix<T>& copy = matrix.Clone();
    copy -= value;

    return copy;
}

template <typename T>
Matrix<T> operator*(Matrix<T>& matrix, T value)
{
    Matrix<T>& copy = matrix.Clone();
    copy *= value;

    return copy;
}

template <typename T> 
Matrix<T> operator/(Matrix<T>& matrix, T value)
{
    Matrix<T>& copy = matrix.Clone();
    copy /= value;

    return copy;
}

template <typename T>
Matrix<T> operator+(T value, Matrix<T>& matrix)
{
    return matrix + value;
}

template <typename T>
Matrix<T> operator*(T value, Matrix<T>& matrix)
{
    return matrix * value;
}



template <typename T>
std::ostream& operator<<(std::ostream& output, const Matrix<T>& matrix)
{
    int columns = matrix.GetColumnsCount();
    int rows = matrix.GetRowsCount();

    output<<endl<<"Matrix with "<<columns<<" columns and "<<rows<<" rows."<<endl;

    int column = 0;
    for(Matrix<T>::const_iterator i = matrix.begin(); i != matrix.end(); i++)
	{
        ++column;
        output<<setiosflags(ios::fixed) << setprecision(2) << setw(5) << setfill(' ') << (*i);

		if(column%(columns) == 0)
		    output<<endl;
	}

    return output<<endl;
}


template <typename T> 
Matrix<T> MatrixMultiplication(Matrix<T>& A, Matrix<T>& B)
{
	if(A.GetColumnsCount() != B.GetRowsCount() /*|| A.GetRowsCount() != B.GetColumnsCount()*/)
		throw INVALIDSIZE;

	int rows = A.GetRowsCount(),
		columns = B.GetColumnsCount(),
		inners = A.GetColumnsCount();
	
	Matrix<T> result(rows, columns);

	for(int row=0;row<rows;row++)
	{
		for(int column=0;column<columns;column++)
		{
			T sum = 0;
			
			for(int inner=0;inner<inners;inner++)
			{
				sum += A.ElementAt(row, inner) * B.ElementAt(inner, column);
			}

			result.SetElementAt(row, column, sum);
		}
	}

	return result;
};

template <typename T>
Matrix<T> ResolveEquation(Matrix<T>& A, Matrix<T>& Y)
{
	Matrix<T> inversed = A.Inverse();
	cout<<endl<<inversed;
	cout<<Y;
	return inversed * Y;
}